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Abstract. We study Lagrangian Perturbation Theory (LPT) and its regularization in the 
Effective Field Theory (EFT) approach. We evaluate the LPT displacement with the same 
phases as a corresponding iV-body simulation, which allows us to compare perturbation theory 
to the non-linear simulation with significantly reduced cosmic variance, and provides a more 
stringent test than simply comparing power spectra. We reliably detect a non-vanishing 
leading order EFT coefficient and a stochastic displacement term, uncorrelated with the 
LPT terms. This stochastic term is expected in the EFT framework, and, to the best of our 
understanding, is not an artifact of numerical errors or transients in our simulations. This term 
constitutes a limit to the accuracy of perturbative descriptions of the displacement field and 
its phases, corresponding to a 1% error on the non-linear power spectrum at k = 0.2/iMpc -1 at 
z = 0. Predicting the displacement power spectrum to higher accuracy or larger wavenumbers 
thus requires a model for the stochastic displacement. 
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1 Introduction 

The detection of the acoustic peaks in the temperature and polarization power spectrum of 
the Cosmic Microwave Background (CMB) anisotropies imply that the seeds for structure 
formation were put in place before the hot big bang, perhaps during a period of inflation. So 
far these primordial seeds are our only fossil record from this early period. 
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The measurements of the CMB anisotropies recently released by the Planck collaboration 
[1, 2] have extracted most of the information available in the CMB about the primordial seeds. 
Further improvements will have to come from measurements of structure at later times in the 
history of the Universe. Gravitational instability makes the primordial fluctuations grow with 
time, making them easier to detect at late epochs. However, the evolution is no longer linear 
and thus on small scales much of the information about the initial conditions is, for practical 
purposes, lost. 

The CMB measurements have produced superb constraints on departures from Gaus- 
sianity of the initial seeds. Although very impressive these constraints fail to reach the clear 
theoretical targets (e.g. [3]), and measurements of LSS are a natural candidate to try to get 
closer to the theoretical thresholds. This has sparked renewed interest in modeling accurately 
the effects of the non-linear evolution of structure in the late Universe, as well as new ideas 
for potential surveys that could reach the required accuracy [4], Accurate modeling of the 
development of structure is also required to extract information about the dark components 
of the matter budget (such as neutrinos) from LSS observations. 

One of the avenues to study the growth of structure is using perturbative techniques. 
Perturbation theory for LSS dates back to the very early days of modern Cosmology e.g. 
[5, 6]. It is extremely successful at calculating correlators at the lowest order or tree level (for 
a complete review of perturbation theory results see [7]). On the other hand, results for the 
first nontrivial correction to tree level results, the “loop corrections”, are less satisfactory. The 
reason for the failure at the loop level is clear. Perturbation theory cannot be used to describe 
the small scales because the series simply does not converge in that regime. Unfortunately 
the errors in the small scales pollute large scale results. This behaviour becomes worse, as 
higher loops are considered, which involve integrals over the spurious perturbative predictions 
on small scales. 

Recent years have seen a quick development of perturbation theory with the introduction 
of the so-called Effective Field Theory (EFT) of Large Scale Structure for both Eulerian and 
Lagrangian descriptions [8-10]. The EFT framework explicitly keeps track of the effects 
of the small scales by keeping explicit counterterms that both fix the mistakes introduced 
by perturbation theory and replace them by the correct level of “leakage” from the small 
scale dynamics. Although the amplitude of these corrections is not known a-priori, their 
functional form is. Thus the EFT of large scale structure generalizes the standard perturbative 
calculation adding new parameters that need to be fitted to simulations or observations. 

In recent years the EFT approach has been used to compute the power spectrum at one 
loop and two loops in a variety of cosmologies, as well as the bispectrum [9, 11-15]. So far, 
the free parameters inherent to the EFT were fitted to simulations, trying to match the power 
spectrum or bispectrum as accurately as possible. This approach is potentially problematic 
as the EFT works best on very large scales where the non-linear corrections are small and 
the sample variance in the simulations is large. So in practice the fitting of these parameters 
was done in the mildly non-linear regime, the same regime were the scheme was being tested, 
allowing the possibility of overfitting. 

In this paper we introduce a novel approach, where we use the EFT to predict the 
displacement field in simulations on a realization by realization basis. Thus we predict the 
actual held in the simulation and not just its power spectrum. By doing so we sidestep the 
cosmic variance issue and directly measure the EFT parameters on very large scales. We can 
then see how these new terms affect the mildly non-linear regime without having to worry 
about overfitting. This approach is closely related to [16], where transfer functions were 
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introduced to quantify the relation between the perturbative predictions and the results of 
simulations. The EFT coefficients are just describing a specific parametrized form for these 
transfer functions. We will generalize [16] by going to higher order in perturbation theory 
and including terms with different spatial structure which are predicted by the EFT. 

This paper is organized as follows: in Sec. 2 we review the Lagrangian perturbation 
theory to establish our notation. We present our simulation suite in Sec. 3, and compare the 
one-loop and higher-loop predictions with simulations in Sec. 3 and Sec. 4. We present the 
evidence for a stochastic displacement in Sec. 5, which constitutes a limit to the accuracy of 
perturbative treatments of the displacement field. We summarize our results in Sec. 6. The 
extensive Appendix gives more details on systematic effects in the simulations, the derivation 
of recursion relations for the perturbative displacement fields, cosmic variance, higher order 
counter terms and transient effects arising from simulation initial conditions. 


2 Review of LPT and EFT: towards a well-defined perturbation theory 


( 2 . 1 ) 


2.1 LPT expansion 

The motion of dark matter particles is described by the Euler-Poisson system: 

dv „ ^ 
at = - Vr$ 

A r d> = 47 rQp. 

Here, v is the physical velocity v = ^j, r is the physical position (r = ax, with x comoving 
position) and t the cosmic time (df = odr, with r conformal time). In the Lagrangian 
picture, we follow the motion of a dark matter particle initially at position q and introduce 
the displacement xp such that x = q + ip. Subtracting the Hubble flow leads to: 


Ip + Uip = 


3 n 3 n 

1 

A cr,(p = 



[Det (Sij + ipij) \ 


( 2 . 2 ) 


Here ‘" ’ denotes the derivative with respect to conformal time, while ‘ * ’ denotes the derivative 
with respect to the comoving Lagrangian coordinate q l , and Ti is the conformal Hubble 
parameter. The above equation breaks at shell crossing, when the determinant vanishes. The 
conformal Hubble parameter is determined by the Friedmann equation: 


V 2 = a 2 Ul 


o° 


(2.3) 


We proceed with a Helmholtz decomposition of the displacement ip = V0 + V Xu. We 
shall call <p the scalar displacement or displacement potential, and u the vector or transverse 
component of the displacement. We make an ansatz for perturbative solutions of the form 

(p = Yju> l ^ and u = Yjr ^ UJ< ' n ’ > ' with: 
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Here <5o is the linear density field, D is the linear growth factor and L n and T n are the LPT 
kernels, which can be computed order by order as we will show shortly [17-20]. Such an 
expansion is exact in the case of an Einstein-de Sitter (EdS) universe, with D reducing to the 
linear growth factor. It constitutes a good approximation for a ACDM universe [7, 17]. 

The goal of this study is to compare the predictions of LPT and EFT for the displacement 
to the non-linear displacement from IV-body simulations. However, the EFT expansion is only 
valid for low values of k, i.e. on large scales, where the cosmic variance in the simulation is 
largest. This cosmic variance introduces relative errors of up to 10% on the power spectrum 
at the largest scales of our simulations (k ~ 5 x 10 -3 h/Mpc), which would make it almost 
impossible to measure the EFT parameters from the power spectrum (see App. C). A simple 
way to reduce this scatter to less than 0.1% is to consider ratios of power spectra (dividing the 
non-linear power spectrum by the linear power spectrum measured from the same realization). 
But an even better way to proceed is to compute the LPT displacement directly on the 
simulation grid, with the same exact initial condition as the one used in the N-body simulation. 
This completely suppresses the cosmic variance. Computing the LPT displacement on the 
simulation grid with the same initial conditions as the simulation will also allow us to perform 
a more stringent test of the EFT, by comparing it to the simulation at the level of the 
displacement held itself, and not only at the level of the power spectrum. 

Since computing the convolutions of (2.4) on the simulation grid would be too compu¬ 
tationally expensive, we will use a real-space formulation of the LPT expansion that is easier 
to evaluate. From the equations of motion presented above, one can solve for the rath order 
displacement recursively (see App. B for the full derivation, following [17, 21]). For the 
scalar displacement, this yields: 
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(2.5) 


where we have introduced 


£-” = )E[d?d:’-d?d:> 

M l,m,n = I At) Am) ,(n) 

OJ ^1*2*3 J1J2J3 


( 2 . 6 ) 


and Ej,j 2i3 is the Levi-Civita tensor. Notice how the scalar part of the displacement on the 
l.h.s. of (2.5) is in principle sourced by both the scalar and vector displacements on the r.h.s., 
through the C and Ai terms. 
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For the vector part of the displacement, we find: 


V x -0 (3) = -n 1 ’ 2 
3 

(2.7) 

n —1 Q 

v x = V 72 — p n p ’ n ~ p 

^ 2n 
p= 1 


with: 

n m,n = V X V (^ W ) ■ (2.8) 

l 

Here again, the vector part of the displacement on the l.h.s. of (2.7) is sourced by both the 
scalar and vector displacements on the r.h.s.. 

Eqs. (2.5) and (2.7) allow to compute recursively the scalar and vector parts of the 
displacement on the simulation grid to arbitrary order in LPT. Notice that the scalar and 
vector displacements are not decoupled: in particular, the vector displacement at third order 
is sourced by the scalar displacements at first and second order, and contributes in return to 
the scalar displacement starting at fourth order. 

From Eqs. (2.5) and (2.7), one can recursively derive analytical expressions for L n and 
T n . Focusing on the equations for the scalar displacement (2.5), we see that there are two 
fundamental vertices corresponding to L and A4. In Fourier space, these correspond to kernels 
which we shall denote by k 2 and ^3 respectively: 


K2(P1,P2) = 1 - 


/ Pi P2 

V P1P2 


2 


(2.9) 


and 


K3(Pl,P2,P3) = 


Det [p \, P 2 , P'sY 
2 2 2 
P1P2P3 


= 1 - 


Pi p 2 

P 1 P 2 


P 2 P 3 

P 2 P 3 


P3 Pl 

P3P1 


+ 2 


(Pi -P 2 )(P 2 -P3)(P3 ’Pi) 


pIpIpI 


( 2 -ld) 

They are the building blocks of higher order vertices. They obey Ko{p, ~p) = 0, Ks(k,p, —p) = 
0. For fixed external momentum k and large loop momentum p —> 00 they scale as 


. . p/k— >00 k 2 

n 2 (p,k-p) - y -3 

pZ 

«3(Pl,P2,fe “Pi -P 2 ) Pl/k ^ 0 °> ^2K 2 {k,p 2 ) 

Pi 


( 2 . 11 ) 


From Eq. (2.5), we derive the LPT kernels L n for the scalar potential in terms of these 
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fundamental vertices K 2 and k 3 : 


L\{p\) = -4 

Pi 


L 2 (.Pl,P2) = ^-2-K2(pi,P2) 


3 

7p?2 


L fa ym) (P\- l P2-,P:i) = -^4“^«3(Pl,P2,P3) 

'V’l23 

L < 3 ™ ym) (p 1 ,p 2 ,p 3 ) = —3— K2(Pl,P2)K2(Pl2,P3) 


( 2 . 12 ) 


7p 


123 


where = p ? ;i + • • • + /L„- 

Similarly, from Eq. (2.7), we obtain the LPT kernels T n for the vector component of the 
displacement: 


T\(.P\) = T 2 {pi,p 2 ) = 0 

r 3 (pi,P2,P3) = -4— Ti(pi)L 2 (P2,P3) (pi ' P23) (pi X P23) 

7*123 


(2.13) 


Eventually, we define the power spectrum Pjj(k) of two random fields r/>W and as: 

= (2ir) 3 5 {D \k + k')Pij(k) (2.14) 


The power spectrum P^ between ith and jth order LPT displacements can be computed 
in terms of the linear density power spectrum Pn n , defined analogously. The 1-loop and 
2-loop power spectra for the the scalar displacement (j> are represented diagrammatically in 
Fig. 1, computed in App. B, and shown in Fig. 2. In what follows, unless otherwise indicated, 
we focus on the scalar displacement and call “ 11 LPT” the scalar displacement up to order n 
(^nLPT = + ••• + </>(")). 

The reader might be more familiar with density power spectra than displacement power 
spectra. For this reason we will often show k^P^Qz) instead of P^(A;), since the former is equal 
to the density power spectrum to lowest order (<5^ = ik ■ r/A 1 ) = —k 2 ^ 1 ^). 

Note that the LPT arising from the kernels in Eq. (2.12) are IR safe by themselves, i.e., 
there is no cancellation of IR contributions as it is the case in Eulerian perturbation theory 
[13]. Thus no IR safe integrand needs to be employed, and the respective integrals can be 
evaluated independently. However, as we will show in detail, the loops are UV sensitive, 
meaning that their value depends on the UV cutoff used, even when this cutoff is already in 
the non-linear regime. This cutoff dependence is unphysical, and makes the LPT predictions 
somewhat ill-defined, as we will show in the following sections. 
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Figure 1. Diagrams for the one-loop (first row) and two-loop (second and third row) contributions 
to the power spectrum of the scalar displacement (j>. The empty squares symbolize the LPT kernels 
L n , while the linear matter power spectrum is represented by the black dots. The explicit formulae 
for the expressions are given in App. B. 



k [TiMpc *] 


k [TiMpc x ] 


Figure 2. Power spectra of the LPT displacements up to 5th order, for one particular realization of 
the linear density field <5 q with cutoff fc m ax = 0.6/iMpc^ 1 . 


2.2 The EFT counterterms 

In the EFT approach, one acknowledges that the mistakes perturbation theory makes on 
small scales will affect the dynamics of the particles on large scales. This can be expressed as 
follows: 

i’ + 'H'tp = -V^^lpt + a ct + a stoc h- (2.15) 

This equation is identical to its LPT counterpart Eq. (2.2), except for a c t, the sum of the 
so-called counterterms, and a stoc h, the stochastic acceleration. 

The counterterms in a c t model the effect of the unknown small scale dynamics on the 
large scales that are being tracked by perturbation theory. Once a few normalization coeffi¬ 
cients (the EFT parameters) are fixed, the Fourier amplitude and phase of a c t are completely 
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determined by the perturbative solution (i.e. ultimately by the linear displacement field) on a 
realization by realization basis. The terms in a c t have the correct spatial structure to absorb 
potential divergencies encountered in the loop calculations. 

On the other hand a stoc h is a stochastic contribution due to the small-scale detail of the 
particular realization considered. Its exact Fourier amplitude and phase cannot be predicted, 
but at the order considered in this paper, only the power spectrum of will be relevant, 

and it will be regarded as realization-independent. 

In [10] the form of the counterterms in the r.h.s. of Eq (2.15) was motivated by con¬ 
sidering the theory for the displacements once smoothed on a sufficiently large scale so that 
perturbation theory is valid. The equation of motion in the EFT can be interpreted as the 
equation of motion of the center of mass of a set of point-like particles, each following the 
LPT equation of motion (2.2). Regardless of the physical origin of these terms, the response 
terms are expanded in powers of the displacement and its derivatives, and to lowest order, 
the resulting expression for the scalar displacement is very simple [ 10 ]: 

^1—loop EFT = (l + otk'' 2 ) (f>^ + 4>^ + 4>^ + 0stoch) (2.16) 

This equation differs from its LPT analog only by the addition of the EFT counterterm 
ak 2 (j>^ (where a is some unknown function of time), and the stochastic term </> s toch- In 
particular, the corresponding 1 -loop power spectrum is given by: 

-Pl-loop EFT = -Pl-loop LPT + 2 ak 2 Pn (2.17) 

Indeed, the power spectrum of the stochastic term is expected to be subdominant compared 
to the EFT term ak 2 P\\ , as shown in Fig. 3 and explained in more details in the next section. 
As we shall show shortly, the EFT counterterm has the right scale-dependence to correct the 
mistake on Pi 3 due to high-fc modes, shown in Eq. (2.24). This way, the EFT term can absorb 
the UV mistake in the LPT power spectrum. Furthermore, the coefficient a will compensate 
the cutoff fc max in the loop integrals, resulting in a cutoff-independent power spectrum. 

Note that we can also understand a as a regulator of loop corrections to the displacement 
field itself (rather than loop corrections to the displacement power spectrum). Namely we 
can close one loop in the third order displacement field (i.e. contract two of the three linear 
density fields) to yield a modified linear displacement field 

0 (3) (fc) D 6 0 (k) I ^L 3 (p,-p,k)PUp) = ( 2 - 18 ) 

The loop is the same one as in P 13 , i.e. it scales as k 2 for large loop momenta. It can 
be regularized by ak 2 . This also makes clear that the regularization of a n-th order field 
(here q in 3 )) requires the leading order counterterm to be a n — 2 -th order field (here <f>^), or 
equivalently that the lowest counterterm counts as a second order field in the power counting. 

In what follows, we call “EFT” or “3EFT” the model cj) 3EFT = (1 + ak 2 )(f)^ + q 
and we shall compare it to the “nLPT” models in terms of agreement with the simulation. 

2.3 Estimates for the sizes of non-linear corrections: EdS scalings vs. ACDM 

When performing a perturbative calculation, it is useful to know how many orders are required 
to reach a desired accuracy level. In the case of the EFT this is also necessary in order 
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Figure 3. Left panel: Power index for each contribution to the power spectrum of </>, as a function of 
the power index of the linear density power spectrum Pi; n . Right panel: Power index of the no-wiggle 
power spectrum from Eisenstein & Hu [22]. 


to understand how many counterterms have to be kept. In this section we summarize the 
standard estimates based on an Einstein-de Sitter (EdS) universe with power law initial power 
spectrum and contrast them with what one gets in ACDM. Since the power law approximation 
is not particularly accurate in the case of a ACDM universe, we will discuss in more detail 
the situation at one-loop order in ACDM in the next subsection. 

We start by reviewing the simple case of an EdS universe, where the initial density 
power spectrum is a simple power law: Pu n (k) = 2TT 2 k~ 3 (k/k n \) n . In this case, each extra 
loop contributes a factor f din q £2 P im (q) (k/k n \) n+3 (when the loop integrals become 
divergent, the situation is slightly more complicated and our discussion here applies to the 
answer after renormalization). One can therefore estimate the size of the 1-loop corrections 
as ~ b(k/k n \) n+3 Pu, and the size of the 2 -loop corrections as b[(k/k nl r+ 3 y P\ 1 , where b 
is some unknown factor. Fig. 3 shows the power index of the various LPT terms in the EdS 
universe, as a function of the power spectrum power index n. 

We compare these EdS scalings to the actual ACDM LPT terms in Fig. 4, where we 
chose n = —1.46 and k n \ = 0.21/i/Mpc to match the amplitude and the slope of the linear 
power spectrum at ko = 0.1/i/Mpc, and we let the factor b float between 0.1 and 10 (colored 
bands). Fig. 4 shows that this approach is too naive. In particular, various LPT terms 
corresponding to the same loop order may differ in size by factors of hundreds on large scales, 
so much so that a 2 -loop term like P 15 can be larger than a 1 -loop term like ^ 22 - It is also 
clear that as we approach k ~ 0.2 — 0.3h/Mpc all terms have roughly the same value. Thus 
one might suspect that at those scales the perturbative expansion fails. We explain this in 
more detail in the case of the 1 -loop power spectrum, in the next subsection. 
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Figure 4. Comparison between the scalings from EdS (shaded areas) and the actual loop contribu¬ 
tions to the power spectrum (solid lines). The scalings correspond to n = —1.4, fc n] = 1-05 h/Mpc 
and b floating between 0.1 and 10. These scalings approximately match the 1-loop terms, but fail at 
reproducing the higher order terms. 


2.4 Size of the one-loop contributions: cutoff-dependence and UV mistake 

We wish to better understand the size of the 1-loop LPT contributions Pi 3 and P 22 , and in 
particular to explain why Pi 3 3 > P 22 • To do so, we decompose and approximate the loop 
integral as 


poo pk poo pk poo 

/ d pf(k,p)= / d pf(k,p)+ / d pf(k,p)~ / d p f(k,p k) + / dp f(k,p^> k), 

■JO JO Jk Jo Jk 

(2 ' 19) 

where in the last step we replace the integrand by its limit when p <C k and p 2 > k respectively. 
This procedure yields the following approximations to Pi 3 and P 22 : 


2P13O) - ^Pn(k) [ e < 5 ,<fc + fe 2 e b,>fc] 

P 22 (k) ~ [2Pn(fc)e* j<fc + rj(k)] , 


( 2 . 20 ) 


where we defined 


e 5,<k 


r](k) 



^P 2p iin{p) OC 
^2 P M(P) OC j 

dp P\m(p) 2 

2 vr 2 p 2 * 


f k dhip 3+ra 

Jo ^ p 

dln P v 1+n 
k 2 vr 2 P 
f 00 dlnp an-r 

Jk 2 tt 2 


( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 


The term eg,<k is the variance of the large-scale density modes with p < k, i.e. it encodes 
tides. The term €,/,>& is the variance of the small scale displacement modes with p > k. The 
contributions of these terms to P 22 and Pi 3 are shown in Fig. 5. Note that the integrand for 77 
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Figure 5. Left panel: contributions to Pi3 from the large-scale density modes (red line) and small- 
scale displacement modes (blue line). Right panel: contributions to P22 from the large-scale density 
modes (blue line) and the p term. The dominant contribution comes from the small-scale displacement 
modes, which only enter the expression of Pi3 and not P 2 2- 




Figure 6. Effect of the cutoff on the loop integral for Pi3 (left) and P 22 (right, same scale). The 
cutoff-dependence makes Pi 3 ill-defined beyond ~ 10 % at k = 0 . 1 /i/Mpc, whereas P22 is very much 
cutoff-independent. 


peaks where the slope of the power spectrum is n = 1/2, i.e., before the wavenumber of matter 
radiation equality at k ~ 0 . 006 h/Mpc. This integral is thus very insensitive to the UV cutoff, 
which is thus the case for P22 as well. It appears that the dominant term for k < 0 . 4 /i/Mpc is 
e^p>k, which is absent in P'22■ This explains why Pi3 P22 for k < 0 . 4 /i/Mpc. Besides, the 

common contribution from es,<k to Pi3 and P22 is multiplied by a larger factor in the case of 
P13. This is why P13 > P22 for k > 0 . 4 /i/Mpc. 

Fig. 5 shows that the low-A: values of P13 are dominated by the contribution 

from high momenta p in the loop. These high momenta cannot be correctly described by 
perturbation theory, and thus bring a mistake in the value of P13. These high momenta in 
the loop also make the result sensitive to the UV cutoff in the integral. On the other hand, 
P22 is very convergent as long as A > k eq , since p(k) is rather insensitive to the slope in the 
UV. Thus the cutoff-dependence of the 1 -loop power spectrum mostly comes from Pi3 (see 
Fig. 6.). This cutoff-dependence is not physical, and means that the LPT power spectrum at 
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(2.24) 


1-loop is at most defined up to ~ 10% at k = O.lh/Mpc. This cutoff-dependence is even worse 
for some higher-loop terms. In particular, Pi 5 calculated with a cutoff of A = 0.6/iMpc _1 
deviates from the calculation without a cutoff by 50% at k = O.lhMpc -1 . The contribution 
to P 33 that arises from the P 33 -n diagram in Fig. 1 inherits the cutoff dependence from P\ 3 . 
The P 33-1 diagram has negligible corrections at k = 0.1/iMpc ^ 1 and P 24 has 10% correction 
at this wavenumber. 

Eq. 2.20 allows us to estimate the scale-dependence of the UV mistake introduced by 
the high momenta in the loop: 

Pv^(k) oc k 2 P u (k) 

P^(k) oc const. 

The EFT counterterms presented in Sec. 2.2 have the correct form to cancel these contri¬ 
butions, and to provide well-defined cutoff-independent predictions. In particular, the term 
ak 2 P\i corrects the cutoff-dependence of P\ 3 , while the stochastic term P s toch corrects the 
(smaller) cutoff-dependence of P 22 . 

2.5 Next to leading order EFT counterterms 

So far we only discussed the leading order counterterm in the EFT. Let us now consider the 
higher order counterterms, that contribute to the two loop calculation. The second order 
counterterm can be derived from a stress tensor T t j. which is constructed from all possible 
second order terms with at least two derivatives of the gravitational potential and two free 
indices [14]: 

Pj = ci 8 \f\d 2 (/>^) 2 + c 2 d i d j (j) 1 d 2 (j)( 1) + czdidi^didj^ + c 4 ^ K) <9;<9 m 0 (1) . (2.25) 

From this we can derive the ansatz for the acceleration counterterm, where we realize that 
only three out of the four terms above are linearly independent. Another derivative yields 
the counterterm for the displacement divergence 

Vo, • a ct = ad 2 (d 2 4 > < ~ 1 ^) + ei<9 2 (d 2 c ^) 2 + e 2 d 2 (d i dj(j ) l ' 1 ' > ) 2 + e 3 di(didj(j)^djd 2 ^). (2.26) 

Reordering the second order counterterms, we define the three linearly-independent operators 


p2,i(fci, k 2 ) —1, 

E 2 ,2{ki,k 2 ) =K 2 (k 1 ,k 2 ), 

E 2 $(ki,k 2 ) = 7 , u 2 y 2 k 2 (Aq, k 2 ), 


(2.27) 


(fci + k 2 

which define the fourth order counterterms by 


k ) = 


d 3 p 
(27r) a 


E 2 ,i{p, k - p)8o(p)5 0 (k - p). 


(2.28) 


These objects are in principle introduced at the level of the equation of motion and need to 
be integrated with the Green’s function to generate a displacement held. These integrals only 
change the unknown prefactors and can thus be reabsorbed into new constants at the level 
of the displacement held. 

Note that two of the three operators, namely E 2)2 and E 2 ) 3 , are also automatically generated 
by solving the equation of motion in the presence of the leading order counterterm ak 2 (j)( 1 \ as 
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we explain in more details in App. E. In particular, these next-to-leading order counterterms 
are generated by either entering the leading order counterterm into /l 1 ’ 1 , i.e. 


k 2(fci, k 2 ) 1 —'^ 2 So(ki)5o(k 2 ) (2.29) 

(« i + «2r 

or by replacing the linear potential in k 4 (j)^ by the second order potential, i.e. 

K2(ki,ki)6o(ki)5o(k2) (2.30) 


Note that these two terms are generated in a fixed proportion, determined by the time de¬ 
pendence of the leading order counterterm a. Note also that E\ is not generated. We will 
see later that all three terms £ 2 , 1 , £ 2,2 and £ 2,3 are needed to cure potential divergencies of 
the displacement field at fourth order. 

There are also new counterterms at third (and higher) order in the fields. We will denote 
them (/>( 3 d) (k) and the corresponding vertices E^j without explicitly deriving their functional 
form. 

We showed above that the leading order counterterm ak 2 cj) can be regarded as a loop 
regulator for the displacement field itself (and not only its power spectrum). Since this 
counterterm is linear in (j), it is a loop regulator for a diagram with a single external leg. The 
diagram for the displacement with one loop and one external leg is obtained by considering 
4>( 3 \ and closing a loop (i.e. contracting two of the three linear density fields). Indeed, we 
found that the ak 2 4> counterterm had the right form to cure the cutoff-dependence of this 
diagram. 

Similarly, the next-to-leading order operators £ 2 , 1 , £ 2.2 and £ 2,3 correspond to coun¬ 
terterms quadratic in </>, and can be seen as regulators for the loops of LPT displacements 
with two external legs (corresponding to terms quadratic in (f >). In order to close a loop, we 
need to add two input fields, which leads us to consider the fourth order displacement: 

1 /* j 3 

^(fci + fco) A -£ 0 (fci)5o(fc2) j (2.31) 

The contributing diagrams, with one loop and two external legs, are shown in Fig. 7 in terms 
of the fundamental LPT vertices k 2 and K 3 . 

Some of them are trivial or reducible counterterms, which are directly related to the 
leading order counterterm that was introduced to cure £13 1 : 


f K 2 (k 1 ,k 2 )K2{ki + k 2 +p, -p)K 2 (k 1 + k 2 ,p) P >fc 24 

PRl oc J - (ftq + k 2 ) 2 - P lin \P) = y^,>^2,2(«l,K2), 


(2.32) 

f K2(fcl,fc2)[«2(fcl +P, -p)«2(fci,p) + {k\ o k 2 )\ p>fc 24 

PP 2 OC J --- (k 1 + k 2 ) 2 - = -^^,>kE 2 ^{k ll k 2 ), 


(2.33) 

Others have a slightly more complex vertex structure, that requires additional counterterms: 


DT f K 3 (k 1 ,k 2 ,p)K 2 (k 1 + k 2 +p,-p) p»fc 12 

pi i°c J - (k 1 + k 2 ) 2 -Pw) = yeV’,>fc^2,2(fci,fc 2 ), 


1 Here we introduce the short-hand notation f 

j p 


f d3p 
J (2n) 3 ' 


(2.34) 
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Figure 7. Diagrams of the reducible and irreducible contributions to the fourth order displacement 
divergence. Note that the vertices in these diagrams are the fundamental LPT vertices n 2 and K 3 
instead of the full LPT vertices Li. The above diagrams arise from writing all possible diagrams for 
</>4 and then correlating intermediate fields, leaving two of the ingoing fields floating. 


PI 2 oc 


[« 3 (fci +P, fc 2 ,-p)«2(fci,p) + Oi k 2 )\ P »fc 12 

- (fci + fc2)2 - = y^ ) >fc^2,3(fci ,fe 2 ), (2-35) 


PIsoc 


n 2 {ki +p, ki - p)K 2 (ki,p)K 2 (k 2 , -p) 


o 


(fet + k 2 y 

48 8 8 

~^zE 2 i(ki, k 2 ) — —E 2 2 (ki,k 2 ) — —E 23 (ki, k 2 ) 

35 35 35 


(2.36) 


We see that all the divergencies can be projected on the basis E 2 ^, which is thus sufficient to 
cure the divergencies. 

Given the leading and next-to-leading counterterms we just derived, we wish to obtain 
the corresponding EFT corrections to the 2-loop power spectrum. At the 2-loop level, the 
LPT power spectra receive corrections that arise from replacing the linear power spectra in 
the one loop expressions by the one loop counterterm k 2 P\[ n , or by replacing the vertices with 
the EFT vertices E, j. as shown in Fig. 8. This leads to the following terms (where the tilde 
denotes a counterterm): 


P 22ti oc. / Pr m (q)Pi in (\k-p\)E 2 (p,k-p)L 2 (-p,-k + p), 
J n 


P~n °c / P hn (p)Pi in (\k-p\)E 2ii (p,k-p)E 2tj (-p,-k + p), 


p i3,i « p \m(k) / P\ in (p)E 3ti (p,-p,k), 

J V 


(2.37) 


P~33,ij « 


Phn(p)P\m{p')Phn(\k ~P~ p'\)E :i Jp, pi ,k - p - p')E 3}j (-p, -pi , -k + p + p') 


p Jp' 


P 33.1 I I Phn{p)Plin(p')Phn(\k - p - p'\)E3,i(p,p' ,k - P - p')L3(-p,-p ,-k + P + p'). 

P J p' 


Note that at the two loop level, only correlators of the second order counterterms with LPT 
terms can contribute: these are the terms P ^ 2 i and P^ i . Thus going to two-loop order in the 
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EFT requires including these counterterms. The auto-power of the second order countert¬ 
erms (i.e. the terms P^ij) an d the auto-power of the third order counterterms (i.e. P 33 ^) 
correspond to at least three-loop diagrams. 

In the tLPT models, which we shall define shortly, any cross spectrum with the LPT 
terms (such as P 22 , and P 13 ?: ) is automatically accounted for by the transfer functions. Thus, 
there is no need to explicitly calculate the counterterms at the two-loop tLPT level. 

In the low-A: regime we have P 22 i oc kr , P 22 2 oc k 6 and P 22 3 oc k 4 , i.e., two of the terms 
scale as P22 while one of them is suppressed by an additional k 2 . For the cubic counterterms 
we expect P 13 oc k 2l P , i = 0,1, 2 . ... These scalings will become relevant later, when we focus 
on the scale dependence of the transfer functions. 

3 Testing LPT and EFT at the field level: moving beyond cosmic variance 
3.1 Simulation Suite: measurements and systematic errors 

We use a suite of Wbody simulations to study the fully non-linear evolution of the cosmic 
density field. The initial conditions for the IVp = 1024 3 particles are set up at ^; n it = 99 
with the publicly available 2LPT code [23]. We subsequently follow their trajectories using 
GADGET-II [24], We use two different box sizes with the same number of particles, the L 
simulation (for “large”) has a box length of 1500/r _1 Mpc and the M simulation (for “medium”) 
has a box length of 500/i _1 Mpc. The cosmological parameters are based on the WMAP7 [25] 
CMB analysis: P m = 0.272, as = 0.81, n s = 0.967. We have 16 independent realizations for 
the L simulation and one for the M simulation. Based on the particle IDs we can reconstruct 
their displacement from the uniform Lagrangian grid 

ip = x — q. (3.1) 

The displacement vector is then assigned to the Lagrangian grid position, the three Cartesian 
components of the displacement grid are transformed to Fourier space and the displacement 
potential is estimated. The number of particles and the box size limit the calculation to 
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wavenumbers smaller than the Nyquist wavenumber k^ y = ^N c /L\) OX , where N c is the number 
of grid cells per dimension (1024 for both simulations) and Lbox is the size of the simulation 
box. This corresponds to A;Ny,L = 2.14/iMpc _1 for the L simulation and &Ny,M = 6.4/iMpW 1 
for the M simulation. 

We calculate the LPT displacements for the same initial conditions, i.e. the phases of 
the perturbative displacements agree with the ones that seeded the fully non-linear IV-body 
simulation. This allows us to perform cross-correlations on a mode by mode basis, rather 
than averaging the perturbation theory and the simulations separately and comparing only 
their auto-power spectrum. This is clearly a much more stringent test of any perturbative 
treatment than just comparing the fairly smooth scale dependence of a power spectrum. 
Matching the power spectrum means that one tries to get 0(1O 2 ) numbers right, whereas 
we are trying to match a O(10 9 ) dimensional vector. In practice, we measure the auto and 
cross-correlations of the various displacements 0 n i, (j)^\ ^ 2 \ and (f>( 5 \ 

In order to compute the LPT displacements on the grid, we use the real-space formu¬ 
lation (2.5) and (2.7). This involves a sequence of products and Fast Fourier Transforms, 
which is computationally more efficient than the convolutions (2.4). In particular, inverse 
Laplacians and derivatives are calculated in Fourier space, but convolutions are computed by 
transforming to real space and multiplying. 

Because this is a non-linear calculation, aliasing is an issue. Coupling n modes with 
wavenumbers k < k cu t 0 fi , we can generate modes up to n fc C utoff- If this wavenumber exceeds 
the Nyquist wavenumber, it folds back to 2k^ y — n ^cutoff and can lead to spurious effects. We 
avoid this by using a cutoff fc cu t 0 ff hr fc-space for our LPT calculation, implemented by setting 
to zero all the modes of the initial linear density held with k > k cu t 0 g. If we want the modes 
k ^ fc cu toff of to be unaffected by aliasing, we need to set fc cu t 0 ff ^ 2/(1 + n)k^ y . 


3.2 Fitting for a 

As a first step in our test of the Lagrangian EFT, we measure the free EFT parameter a that 
appears at 1-loop order, as defined in Eq. (2.16). We fit for the coefficient a in two ways: by 
minimizing the variance between the non-linear held and the theoretical model on the one 
hand, encoded by the spectrum -P e rror of </ error = </> nl — •/’modeb and by adjusting the model 
power spectrum to the non-linear power spectrum P n i on the other hand. The resulting ccfc 
are obtained as follows, and shown in Fig. 9: 


minimizing P error = (|</ nl - ^ 3LPT - ak 2 c /) (1) \ 2 ) 

< 

minimizing (P n i — P 3 lpt — 2 aA: 2 Pii) 2 


Pi X nl Pi 


QWror — 


lx3LPT 


Q( n l — 


k 2 Pu 

1 Pnl ~ -P 3 LPT 


k 2 P\ 


11 


(3.2) 


In both cases, we first determine the value of ak for each k-bin independently, and 
estimate the error bars on these values by looking at the scatter across the independent 
realizations. Keeping the hrst term ak 2 in the EFT expansion is only supposed to be a valid 
approximation at low k, where the higher order terms are small. We would therefore expect 
the ak curves shown in Fig. 9 to look like a constant at low k, and then to deviate from that 
constant at higher k. 

Instead, the curves we obtain seem to deviate dramatically from a constant at the 
smallest k, for the following reason. At very low k, our measurement of a is extremely sensitive 
to systematic errors in the simulation, since an error of x% on the non-linear displacement 
translates into an error of x/k 2 % on a. The grey domains on Fig. 9 represent the variations 
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in a caused by variations in the non-linear displacement of 0 . 01 %, 0 . 1 % and 1% respectively. 
On large scales (k < 0.05/iMpc _1 ), the dramatic upward or downward trends seen in a (right 
panel of Fig. 9) are the reflection of the ~ 0.05% variations in the non-linear power spectrum 
arising from different IV-body settings discussed in App. A. 

The values of a obtained also deviate from a constant at high k, which is expected due 
to higher order terms in the EFT expansion. Some of these terms - the ones that correlate 
with - affect both methods in the same way. This is the case for the Pi 5 term, which we 
will discuss in more detail below in Sec. 4.2, and the higher order order EFT counterterm of 
the form . A naive estimate for the latter term is given by a 2 * /c 4 Pn, which leads to a 

5% correction to the inferred value of a at k ~ 0.14/iMpc _1 , where a itself leads to a 10% 
correction to the non-linear power spectrum. 

However, some other higher order EFT terms affect only the value of a & obtained from 
the non-linear power spectrum. These are the terms that are not of the form Pi*. Indeed, 
considering only the LPT terms, we expect 


terror — 


Pllx 1 — Pll — P] 


13 


k 2 P 1 


Q( n l — 


Pnl — Pll — P22 — 2 Pl 3 


11 


2k 2 P\ 


11 


and therefore 


P33 + 2P24 

«nl = terror + „ „ - + 


(3.3) 


(3.4) 


2k 2 P u 

These correction terms are suppressed on large scales, but sufficient to explain the percent 
level difference between the a*, obtained from the two methods, shown in Fig. 9. Thus, 
as soon as the a.k from the two methods differ, we know that higher order corrections are 
important, and the estimator from fitting to the non-linear power spectrum should not be 
trusted anymore. At the same time, two-loop corrections to the power spectrum will become 
important. Note that the auto-power spectrum of the stochastic term will only affect the 
value of a n i and not a e rror- 

From Fig. 9, the values of a inferred from both methods agree for k < 0.03h/Mpc, 
and disagree for k > 0.03/i/Mpc. This shows that for our purpose, where we focus on very 
large scales, there is no overfitting issue: choosing a. to get the best displacement held or the 
best power spectrum on these scales is equivalent. However, this means that one should not 
determine the value of a by matching the power spectrum on scales k > 0.03/i/Mpc. 

Picking the value of a around k = 0.03/iMpc -1 , we find a = —2.46 ± 0.2(h - 1 Mpc) 2 , 
where the error bar represents the deviations between simulations run with different settings, 
and therefore corresponds to a systematic error. Remember also that the exact value of a 
is dependent on the cutoff used for the LPT calculation. Here we used fc cut off = 0.61/t/Mpc. 
This can be translated to q:a=oo = —3.0 ± 0.2(h - 1 Mpc ) 2 for /c cu toff = 00, which corresponds 
to a 1% correction to the power spectrum at k = 0.06/iMpc -1 . 2 The shaded areas in Fig. 9 
show that our value of a allows to reproduce the non-linear power spectrum with an accuracy 
of 1% up to k ~ 0.15/i/Mpc. 

The same measurement at z = 0.5 yields the value a = —0.73 ± 0.07(/i _1 Mpc) 2 (see 
Fig. 10), corresponding to a \ =00 = —1.13 ± 0.07(h _1 Mpc) 2 for fc cu toff = 00 . However, at 
higher redshift, the absolute value of a drops below the systematic error in the simulations, 
and we are no longer able to measure it. 


2 This k 2 coefficient is not captured by two loop corrections. In fact there is a positive fc 2 Pn contribution 

from Pis, which would thus require an even more negative value of a after the two loop power spectrum has 

been included. We will come back to this issue in Sec. 4.2. 
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Figure 9. Left panel: Values of ctk obtained from fitting to the power spectrum and from minimizing 
the error power spectrum for a cutoff fc cu toff = O.bl/iMpc -1 , extracted from the same realization. The 
two methods agree on large scales but start to differ at percent level for k > 0.03/iMpc^ 1 . Right 
panel: Values of at from minimizing the error power spectrum, for the different simulation runs and 
^cutoff = O.Gl/iMpc -1 . This shows that the systematic error in our measurement is sizeable for the 
smallest k. The shaded areas indicate the domains in which changing a does not affect the non-linear 
power spectrum by more than 0.01%, 0.1% and 1% respectively. 



Figure 10. Values of a*, from minimizing the error power spectrum at redshift z = 0.5, for the 
different simulation runs. The shaded areas indicate the domains in which changing a does not affect 
the non-linear power spectrum by more than 0.01%, 0.1% and 1% respectively. 

Cutoff-dependence Loops are integrals over all momenta, and therefore introduce a mis¬ 
take due to the fact that the high momenta cannot be described correctly by perturbation 
theory. The role of the EFT coefficient is to compensate this mistake, and thus the exact 
value of a depends on the high-A: cutoff used in the perturbation theory: 

«(Al) = <*(A 2 ) + | . (3.5) 

This cutoff-dependence is shown in Fig. 11. In particular, sending the cutoff to infinity in 
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Figure 11. Left panel: Cutoff-dependence of the EFT coefficient a. The red data point is measured 
from the L simulations at z = 0. The blue band corresponds to the fitted value plus the correction 
from Eq. (3.5). The width of the line is given by the errorbar on the data point. Right panel: Measured 
value of a, as a function of redshift. The data points have been extrapolated to fc C utoff = oo. The blue 
line corresponds to the naive scaling oc 1/ k^\ which we would expect for a scaling universe. Here k n \ 
is defined by k^P(k n i, z)/(2 tt 2 ) = 1. 


the loop integral results in a non-zero value a = —3.0 ± 0.2 (Mpc /h) 2 . The coefficient a can 
be set to zero by choosing /c C utoff — 0.15/i/Mpc, but this has nothing physical, and does not 
imply that the EFT counterterm was not needed. 

Redshift-dependence For an EdS universe with power law initial power spectrum of index 

4 

n, the relation between length and time scalings is = A ” +3 [12]. We expect the EFT 

„ 8 4 

coefficient a to scale as its dimension, i.e. length , i.e. as r n + 3 oc a n + 3 . Measurements of 
aA=oo a t different redshifts are shown in Fig. 11. The time-dependence of qa=oc appears 

9 9 4 

compatible with that of length' oc l/k^{a). and with D n + 3 with n = — 2 . 

3.3 One-loop power spectrum 

We compare the non-linear power spectrum from LPT and EFT to the simulation output. 
Fig. 12 shows that the 1-loop EFT power spectrum matches the simulation to 1% accuracy 
up to k = 0.15/iMpc _1 , compared to k = 0.05/tMpc -1 for the 1-loop LPT power spectrum. 
This is a factor of three improvement in the maximum wave vector giving 1% accuracy. 
Let us stress again that the EFT coefficient a was not fit at the maximum wavenumber of 
validity k = 0.15/iMpc -1 but on larger scales. Notice that the 2-loop LPT power spectrum 
is worse than the 1-loop LPT power spectrum, which is not surprising given that it contains 
UV-sensitive terms that haven’t been corrected by the appropriate EFT counterterms yet. 

It is striking from Fig. 20 how good the agreement is at low k, completely devoid of 
cosmic variance. This shows that the EFT model works realization by realization, instead of 
only reproducing the mean power spectrum. 

Thus, the third order EFT displacement provides a very good fit to the non-linear 
displacement on weakly non-linear scales, where the expansion in powers of (k/k n {) is valid. 
However, one should not evaluate it at high k. In particular, computing the root mean square 
displacement from the EFT model yields wrong values, because this calculation relies on the 
displacement at high k for which the EFT term oc k 2 diverges. 
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Figure 12. Left panel: Relative difference between the non-linear power spectrum from the simulation 
and from LPT (Zel’dovich approximation, 1 and 2-loop) and EFT (1-loop). The 1% error domain is 
the shaded grey band. The maximum wave vector with accuracy of 1% is improved by a factor of 
three from 1-loop LPT to 1-loop EFT, from 0.05 hj Mpc to 0.15 h/ Mpc. The 2-loop LPT worsens 
the agreement to simulation, compared to 1-loop LPT. The shaded magenta region indicates the 
scatter we would get due to cosmic variance without the LPT calculation on the simulation grid: 
this measurement has negligible cosmic variance. Right panel: Power spectrum of the error on the 
displacement field. Adding the second and third order to the first order displacement improves the 
agreement at the level of the displacement field on large scales (k < 0.1/iMpc -1 ). Including the EFT 
counterterm at 1-loop further improves the agreement, by correcting the UV mistake in However, 
going up to fifth order in LPT worsens the agreement, as expected for an asymptotic series, because 
of the UV mistake that is not corrected by EFT counterterms. 


3.4 Relative importance of the various EFT terms 

As Fig. 12 shows, the EFT provides a good fit not only to the non-linear power spectrum, 
but also to the displacement field itself. However, in the case of the EFT power spectrum, 
the contribution from (j)^ (i.e. the term P 22 ) is negligible compared to the contribution from 
(i.e. the term P 13 ). One might therefore wonder about the relative importance of the 
non-linear terms <f>^ 2 \ ak 2 q present in the EFT model: do they contribute equally? 
Is the second order displacement helping at all in the agreement with simulation? 

The answer to these questions can be visualized as follows. The displacement fields 
0 nb are functions of the wave vector k, i.e., they are defined for each of the 

iVmodes modes in our simulation box. They can thus be understood as very high dimensional 
vectors {<j>(ki)) i=1 N . We can then interpret {(j> a \^>b) = {^a^b) as a scalar product and 
(|^ a | 2 ) as the corresponding squared norm on this vector space. Intuitively, with this scalar 
product, two displacement fields are aligned if they are perfectly correlated, and orthogonal if 
they are completely uncorrelated. This allows a graphical representation of the displacement 
fields on the basis (j)^ 2 \ <^ 3 ') of the LPT terms. This basis is not orthogonal (e.g. 

(c( 3 )) 7 ^ 0 ), so we shall instead use the orthonormal basis <^ 2 ^, , deduced 

from </>( 2 \ (j)^) through the Gram-Schmidt orthogonalization process. Fig. 13 shows 

the graphical representation of the EFT terms as well as the non-linear displacement. Fig. 13 
makes it visible that the contribution of f>: 2 > to reducing the error terror is more important 
than that of (j)^\ as one would expect for a well-behaved expansion. It also shows that even 
though is a smaller term than (f>^ (i.e. ||0( 3 )|| = vPi .3 < V^2 2 = ||^ ( ' 2 ^||)) it brings a 
larger contribution to the non-linear power spectrum, because it is more “aligned” with cj>^ 
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Figure 13. Vector space representation of the LPT terms </>0) (blue), c/4 2 ) (green), q i4 3 l (red), q 
(magenta), the EFT term ak 2 (f>^ (cyan) and the non-linear displacement (f> n \ (black), at k = 0.045 
h/Mpc. The shaded grey sphere corresponds to a displacement error such that ||terror|| = l%||<^ni|| 
(i.e. V-Perror = Pn\) ■ This representation shows that </4 2 ) is indeed a bigger term than <(7 3 ), and 
contributes more than to reducing <p errol . It also shows that the effect of <f > < - 2 - ) on the non-linear 
power spectrum is smaller than that of 4 >^ 3 \ because q^ is more “aligned” with </>d) than q 7 2 ) is. 

(i.e. (^W(/>( 3 )) = Pi 3 > 0 = {(j. !>W(/>( 2 ))). In conclusion, the second order displacement q^ is 
a larger term than the third order displacement (f )^, and is crucial to the agreement between 
EFT and simulation at the level of the displacement field. However, because q in 2 ' is orthogonal 
to (i.e. not correlated with) the larger term (j)^\ its contribution to the displacement power 
spectrum is negligible. 

3.5 On the overfitting issue 

In this section, we fitted for the EFT coefficient a, and compared various models of the 
displacement field to simulation, at the level of the field itself, and not only its power spectrum. 
This methods allows us to address the issue of overfitting. 

Indeed, one might a priori be concerned that the EFT expansion for the displacement 
field might be incorrect, but still give the right displacement power spectrum due to the 

adjustable EFT parameters. For example, the (artificial) model </> m odel(&) = 
clearly predicts the correct power spectrum P n \ for all k, but as shown in Fig 14, it does not 
correspond to the true non-linear displacement (j> n \. One might worry that the EFT model 
might be similar, in that it would give an accurate non-linear power spectrum thanks to its free 
parameters, but be a poor description of the displacement field. However we did not choose 
the EFT coefficient a by requiring P mo dei to be close to P n i, but instead by minimizing the 
error power spectrum Herron which is the most stringent requirement. Indeed, a small P ervor (k) 
means that |terror (k)\ 2 = \4> n \(k) — </>model (&)| 2 is small, implying that the Fourier component 
of the scalar displacement is correctly described by the model. In contrast, having P n \ close 
to Hmodel simply means that \(f) n \(k)\ 2 — | pmodei (&)| 2 is small, i.e. the Fourier components 
have the same amplitude, but might have different phases. Therefore, by minimizing -Perron 
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Figure 14. Comparison between the error power spectra for the artificial model </> mo dei(&) = 
\J pT i(k) (&) an d </> 3 LPT- Even though the former gives exactly the right non-linear power spectrum, 
its error power spectrum is much larger, meaning that it is a wrong model for the displacement. 


we are guaranteed to choose the value of a that provides the best displacement field possible, 
which prevents the risk of overfitting. 

In practice though, because our measurement of a is largely limited by systematic er¬ 
rors in our IV-body simulation suite, the difference between the value of a obtained from 
minimizing the error power spectrum and from fitting to the non-linear power spectrum is 
small compared to the final uncertainty on a. This was not a priori obvious, and shows how 
sensitive this measurement is to systematic errors on the largest scales of the simulation. 

4 Transfer functions and ‘tLPT’ 

4.1 Optimal linear model from LPT — including higher-loop terms 

In the previous section, we measured the EFT coefficient a, and compared the EFT model to 
the various nLPT models. We found that the EFT significantly increases the range of validity 
of perturbation theory. But we also wish to understand how close to optimal the EFT model 
is. To do so, we compare it to the “ntLPT” models (for “LPT with transfer functions”), i.e. 
the LPT models for which the LPT displacements are multiplied by free functions of the 
modulus k = |fc| of the wavenumber (the “transfer functions”, see also [16]): 

4>ntLPT{k) = ai(A:)0 (1) (fc) + ... + a n (fc)</> (n) (fc). (4.1) 

For each of these models (LPT, EFT, tLPT), we compute the displacement error terror = 
0nl — ^model• The transfer functions di(k) of the tLPT models are chosen so as to minimize the 
error power spectrum Terror for each /c-bin. Notice that the transfer functions are not chosen 
so as to match the non-linear power spectrum. The tLPT models correspond to a lower bound 
on Terror) and thus allow to assess how close to optimal the LPT and EFT models are. 

Besides, the transfer functions allow to effectively include certain higher order LPT and 
EFT terms without having to compute them explicitly. For example, for the ltLPT model 

0itLPT = the result of minimizing T err0 r yields ai = P p* nl , so that TitLPT = p*” 1 — 
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Figure 15. Transfer function ai of the ltLPT model. It corresponds to a\ = ^ pl ~ 1 + ak 2 + 
which indicates that the ltLPT power spectrum effectively includes the Pl 3 term, even though q 
is not explicitly included in the ltLPT model. 


Pn + Pi 3 + ak 2 Pn + ... . Here, the term P 13 is implicitly included in the ltLPT power 
spectrum, even though evaluating the ltLPT displacement did not require computing 

In turn, the LPT and EFT can predict what these transfer functions should be, as 
shown in Fig. 15. We see that the 3LPT model for ai overpredicts its amplitude due to the 
UV sensitivity of P 13 , whereas the 3EFT prediction, taking into account the counterterm 
associated with a, is accurate to 1% even beyond k = 0.1/iMpc -1 . 

This also generalizes to higher order: as soon as the term atf is included, the LPT 
contributions Pij to the non-linear power spectrum for all j are automatically included. This 
is illustrated in Tab. 1. 



Table 1 . List of contributions to the LPT power spectrum at tree order and various loop orders. The 
columns collect the terms that are effectively included by the transfer function (ai, a 2 , a 3 , 04 and 05 
respectively). For example, this shows that 3tLPT contains all of the 2-loop terms, including Pi 5 and 
P 24 , even though <j >4 and <f >5 are not computed explicitly in </>3tLPT- 


In practice, we fit for the free coefficients a.j(fc) for each k -bin independently, by mini¬ 
mizing Perror(^) = ( terror | 2 ) = (|</>nl - a i4>i | 2 )- If we interpret {(j)a\<t>b) = {K^b) as a scalar 
product and (|^>a| 2 ) as a squared norm as before, we see that minimizing the norm of the 
displacement error terror (be. minimizing P e rTO r) amounts to making the displacement error 
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Figure 16. Left panel: cross-spectra between the LPT terms and the non-linear displacement needed 
to make comparisons at the field level, and to fit for the transfer functions (defined later). These were 
computed for one particular initial condition. The shaded area indicate the uncertainty due to cosmic 
variance. The cosmic variance is visibly much higher for P n i X j with i even, as explained in App. C. 
Right panel: Ratios of the cross spectra from the left panel and the lowest order LPT prediction. The 
cosmic variance on the ratios is significantly reduced. 


orthogonal to the LPT components (f>i (i.e. (f> 

X 2 = <knl-E a i^ W | 2 > minimal ^ = - 


error |#) =0): 


- 2 ( 0 «| 0 error ) = 0 . 

(4.2) 


Since the various <j)W do not form an orthogonal basis (i.e. ^ 0, i.e. (f>^ and 

can be correlated), the value of the transfer function depends on the maximum order in¬ 
cluded in tLPT: for instance, the value of a\ in ltLPT and 3tLPT differs by Pvz/Pu- This can 
be avoided by defining an orthogonal basis through the Gram-Schmidt orthogonalization 
process, applied to the LPT displacements (j>^: 


^ = </>& 





(0( i + 1 )|^W- L ) 

(^CO-L^W-L) 


(fU) 1 - 


(4.3) 


This way, we can define “orthogonal” transfer functions a^~, which are independent of the 
order in tLPT: 

^ntLPT = + ••• + «n^ (n)_L (4.4) 

Note that this does not affect the value of the tLPT displacement: it is a decomposition of the 
same tLPT model in the basis of the instead of the (</>^). However, this expression 

will be useful in estimating the residuals <f n i — ^ n tLPT) as we shall see later. 

As we have seen, fitting for the transfer functions requires knowing the cross-spectra 
(</>W| (f)^)) between the various LPT terms, but also the cross-spectra (4>^\4>nl) between the 
LPT terms and the true non-linear displacement. For reference, we show the latter in Fig. 16. 

The measured orthogonal transfer functions a^, , 0:3 and CI 4 at redshift z = 0 are 

shown in Fig. 17. The non-orthogonal transfer functions of 4tLPT at redshift z = 0 and z = 2 
are shown in Fig. 18. The measured transfer function a n or a^ becomes increasingly sensitive 
to potential systematic uncertainties in the simulation as n increases. 
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Figure 17. Orthogonal transfer functions at redshift 0. The shaded area corresponds to the domains 
where varying a,i changes the contribution to <^4tLPT by less than 0.1%. 




Figure 18. Non-orthogonal transfer functions at redshift z = 0 (left panel) and z = 2 (right panel). 
The dashed lines are from the L simulation and the solid lines from the M simulation. At redshift 
z = 0 the transfer function <34 of the fourth order displacement on large scales is close to zero, but it 
approaches unity in the M simulation at z = 2. Such a behavior can be expected if the k° component 
of P46/P44 is of order unity. Note also that the transfer function a\ on the Zel’dovich displacement 
remains close to unity up to non-linear wavenumbers. 


4.2 Transfer functions in the context of the EFT 


We mentioned above that 3tLPT includes all the terms present in a consistent two-loop EFT 
calculation, including the next-to-leading order counterterms. Here we make the connection 
between EFT and tLPT explicit, by showing that the transfer functions are correctly described 
by LPT power spectra and their EFT counterterms on large scales. For clarity, we focus on 
the orthogonal transfer functions a^, which are independent on the order used in tLPT. 

The EFT predicts the first order orthogonal transfer function as follows: 



Elxnl 


~ 1 + 
~ 1 + 


Pi 3 

P 11 

P\ 3 

P 11 


+ ak 2 
Pi 5 


+ 


E11 


+ a'k 2 . 


(4.5) 
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Figure 19. Left panel: Transfer function aj~ of the ltLPT model. Adding Pi5 modifies the value of 
a, as expected, and improves the agreement. Right panel: Transfer function a2 of the ltLPT model. 
We clearly see a percent level deviation on the largest scales, that is accounted for by adding P24/P22 
to the model for this term. The latter however overpredicts the enhancement in the mildly non-linear 
regime, which is in turn fixed by the EFT counterterms P 2 y. As we pointed out before, they lead to 
k° and k 2 corrections through P^fPii- 


Here, the first line corresponds to the 1 -loop EFT prediction for a^, where a is the EFT 
coefficient measured above. The second line corresponds to the 2 -loop prediction (see Fig. 19 ). 
Since the LPT terms such as P13 and Pi5 are cutoff-dependent and have potentially wrong 
UV contributions, they are associated with a counterterm. Thus a and ol differ by the k 2 
coefficient of P15/P11 which is approximately 3 /i _2 Mpc 2 for a cutoff of A = 0 . 6 hMpc -1 . 
Beyond that, P13/P11 and P15/P11 have k 4 contributions which are cutoff dependent, and 
should be corrected by a counterterm of the form / 3 fc 4 . However, we found that including 
such a counterterm does not improve the prediction for af. and we therefore do not include 
it in the following discussion. 

Similarly, the deviations of the transfer function for the second order displacement can 
be modeled by a combination of the next-to-leading order LPT term and the corresponding 
EFT counterterms 


P24 , P 


a 2 ~ 1 + if" + ~p^' 

-T22 422 


( 4 . 6 ) 


Note that P24/P22 starts as k° (actually predicting deviations of 1 % on large scales at z = 0 ). 
This means that the transfer function on <f>2 can deviate from unity even on very large scales. 
The cutoff dependence of P24 needs to be captured by the corresponding counterterms - 
P22/P22 is the sum of three next-to-leading order counterterms with their free coefficients. 
As discussed in Sec. 2 . 5 , they scale as k° and k 2 for small wavenumbers. We show the EFT 
description of the second order transfer function in Fig. 19 . P24 captures the large scale offset 
quite well, but overpredicts the scale dependence of the transfer function. This mistake can 
be corrected by adding the quadratic counterterms in the form of a k° term and a k 2 term 
with free coefficients. 

More generally, the transfer functions beyond the one on do not necessarily go to 
unity on large scales: they can be renormalized at the k° level by higher order LPT and 
EFT corrections. The reason for the distinction between aq and a>i is that the LPT kernels 
already ensure mass and momentum conservation in the coupling between initial modes and 
can thus be modified at the k° level. For the linear field (f^ l \ corrections need to explicitly 
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Figure 20. Comparison between 1-loop LPT, EFT and tLPT, in terms of their agreement with 
simulation for the power spectrum (left panel), and the displacement field itself (right panel). The 
EFT improves on the LPT, and is close to tLPT, with the latter performing slightly better on the 
power spectrum. This shows that the 1-loop EFT cannot be outperformed by any expansion of the 
form ai(k)(f>^(k) + ci 2 {k)(j)( 2 \k) + a 3 (k) 4 >^ (k) by more than a factor of two in the wavenumber up 
to which we can trust the theory at percent level. 

conserve mass and momentum and thus need to start as k 2 q [6, 26]. 

4.3 Optimality of the EFT at 1-loop 

As explained earlier, we compare the EFT model to the tLPT models to assess how close 
to optimal the EFT model is. In terms of the error power spectrum, which quantifies the 
agreement with simulation at the level of the displacement held, the 1-loop EFT and tLPT 
both have Terror ^ l%T n ] up to k — 0.2hMpc _1 . This shows that the EFT displacement 
is very close to the 3tLPT displacement, showing that the EFT model is close to optimal 
at 1-loop order: no expansion of the form ai (k)t/>^(k) + a 2 (k)(j)^(k) + a 3 (fe)0 (3) (fc) can 
significantly outperform it. In terms of the non-linear power spectrum, the 1-loop EFT 
provides a 1% fit to the power spectrum up to k > O.lhMpc -1 , compared to k ~ 0.2/tMpc -1 
for the 3tLPT model. One might be able to achieve this factor of two extension in the range 
over which the power spectrum can be described at the 1% level using 2-loop EFT. 

5 The stochastic term 

5.1 A floor in the error power spectrum 

Detection of the stochastic term We wish to understand the error power spectra 
for the various tLPT models. We make use of the orthogonal basis defined above. In 
this basis, the transfer functions af{k) are independent of the tLPT order used (i.e. ajq(k) 
is the same for ltLPT, 2tLPT, etc). Since (j> 4tLPT is our best model for the true we will 
write: 

<£nl ~ Oi L 0 (1) " L + a2^ 2)± + «3 ^ (3)± + «4 </ )(4) " L + ^stoch (5.1) 

Thus we can estimate the displacement errors for ltLPT, 2tLPT and 3tLPT as follows: 

terror' 1 ' - a 2 + a 3 + af + <?Wh 

- 03 ^ (3)± + «4 ^ (4)± + </Wh (5-2) 

terror" 1, - a A + 4toch 
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Figure 21. Error power spectrum for the various tLPT models (solid lines), compared to the naive 
expectation of equations (5.1) and (5.3) (dashed lines), relative to the non-linear power spectrum. 
The naive expectation underestimates the measurement, which indicates the presence and the exact 
size of the stochastic term P s toc h- 

For a well-behaved perturbative series, one would assume that the various terms on the 
r.h.s. of (5.1) are ranked in decreasing order. Keeping only the dominant term and neglecting 
the stochastic term then leads to the following estimate for P e rror for the tLPT models: 

« T -[4] 2 ^ 2± 

« T -[4] 2 P3T3T (5-3) 

- [4] 2 lw 

These estimates are shown in Fig. 21 (dashed lines). The fact that they respect a^P^^ 0 > 
a 3 P):) tho > a^P^ ho shows that this expansion is well defined: higher order terms are indeed 
smaller. However, Fig. 21 also shows that these estimates for the are much lower than 

the measured error power spectra (solid lines). There is clearly a floor in the measured error 
power spectra, which we associate with the stochastic displacement </> s tock- Fig. 21 shows that 
the correct ranking of the terms in the 4tLPT model on large scales (k < 0.1/iMpc -1 ) is not 
that of (5.1), but instead 

4>n\ = aj L 0 (1)± + 0,2 </> (2)± + O 3 <j 6 (3)± + </>stoch + of (/> (4)± , (5.4) 

i.e., the stochastic term is not negligible and even exceeds the amplitude of the fourth order 
displacement. These findings indicate that the error power spectra should scale according to: 

^errof T - [o^] 2 P 2 _L 2 ± + ^3 ] 2 P 3 _L3_L + -Pstoch 

■^error* T - [ a 3 "] 2 ^3±3± + f’stoch ( 5 - 5 ) 

p^r ^ ^stoch 

Indeed we can define a unique function of k, the stochastic power spectrum, for which all 
three of the above equations are satisfied. We show the stochastic power spectrum from the 
L simulation as the solid black line in Fig. 22. 
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Figure 22. Ratio of the error power spectrum for ltLPT, 2tLPT and 3tLPT and the non-linear 
power spectrum after the stochastic term has been removed (dashed lines, same as Fig. 21). The 
stochastic term P stoc h (black solid line) corresponds to P error = l%P n i at k ~ 0.2 h/Mpc. Beyond 
3tLPT, the new LPT terms have contributions much smaller than the stochastic term, and therefore 
do not bring any significant improvement on the error power spectrum. 


The term pstoch is by definition not correlated with the LPT terms, and acts as a 
noise that limits the accuracy achievable with tLPT (see Fig. 22): because of this term, the 
improvement in P e rr0 r from 2tLPT to 3tLPT is limited, and there is no improvement at all 
beyond 3tLPT. The ragged features of Pstoch visible in Fig. 22 are likely a consequence of 
systematic errors in our IV-body simulation. They correspond to very small errors on large 
scales: P err0 r ~ 10 _6 P n i at k ~ O.OlhMpW 1 . In Fig. 23, we show the stochastic power 
spectrum for the L and M simulations, and its time dependence in the M simulation. We 
clearly see that the stochastic power spectra agree between the two simulations. 
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Figure 23. Power spectrum of the stochastic term. The red points show the stochastic term in the 
M simulation at redshifts z = 0,0.5,1, 2 from top to bottom. The stochastic term of the L simulation 
at z = 0 is shown by green triangles and agrees perfectly for k < 0.6hMpc -1 except for a small upturn 
on large scales. The thick shaded line shows the systematic error on the L simulation at z = 0, which 
exceeds the stochastic term on the largest scales, where the two simulations disagree. The red lines 
show the phenomenological model for the stochastic term from Eq. (5.10). 


Scale-dependence: On large scales, we can fit the stochastic term in Fig. 23 by 

h « 5.7 X 10* (ft-'Mpc) 3 D w ( , M *p c -l) 4 • 


(5.6) 


On the largest scales, we thus find that P s toch is independent of k. This is what we expect from 
mass and momentum conservation: for a mass and momentum-conserving perturbation, the 
lowest order correction to the density field is expected to scale as k 2 [6, 26], which corresponds 
to a ^-independent correction to (j). The stochastic term predicted by the EFT comes from 
a small-scale reshuffling of the matter, which is not describable in terms of the large-scale 
displacements, but conserves mass and momentum on large scales. The scaling observed in 
(5.6) is thus consistent with it being the stochastic term expected in the EFT framework. 
More specifically, in a scaling Universe we expect 


A 


2 

stoch 


k 3 (k 4 P stoch ) 

2-7T 2 



(5.7) 


and self-similarity dictates k n \ oc a~ 2 ^ n+3 \ where n is the slope of the initial power spectrum. 
The fitted time dependence D 10 is reproduced by a slope of n = —1.6, which corresponds 
to the slope of our input power spectrum at k = 0.1/iMpc -1 . Furthermore, from A 2 toch = 
(k/k n \Y we can deduce k n \ ~ 0.32/rMpc -1 at z = 0. Due to the steep scaling ( k 7 ), an order 
one prefactor does not change k n \ significantly. We show the time dependence of the power 
law part of the stochastic term and the time dependence of the scale where it amounts to a 
1% change in the power spectrum in Fig. 24. 


Amplitude of the stochastic displacement To get an idea of the size of the stochas¬ 
tic term, we use the measured power spectra to infer the corresponding root mean square 
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Figure 24. Upper panel: Time/redshift dependence of the large scale amplitude of the stochastic 
term. Lower panel: Time/redshift dependence of the scale, where the stochastic term amounts to 1% 
of the total power. 

displacement 

2 If d 3 fc . If d 3 k 2 
^.„ch = 3/(^3 <*!*> = 3 / (5.8) 

and find a rms stochastic displacement <T stoc h — 0.8 Mpc/h (compared to a n \ ~ 6.0 Mpc/h for 
the full non-linear rms displacement). The integrand in Eq. (5.8) is peaked at k ~ 0.6 fiMpc -1 , 
i.e. far beyond the range of scales where the stochastic term follows the k 4 scaling. 

We do not expect perturbation theory to be able to capture the motion of particles within 
halos. These motions will thus contribute a (probably significant) fraction of the stochastic 
displacement. A crude estimate for this source of stochastic displacement can be obtained 
by assuming that particles within a halo of mass m have a root-mean-square displacement of 
roughly two virial radii i?. v j r . and to average over all halos using the halo mass function: 

<4och.,h = J dm ^ m ( 2 - R virM') / J ( 5 - 9 ) 

Integrating over all masses, this estimate yields er atoc h. motion in halo — l.l/i _ 1 Mpc, similar 
to crstoch• The integral peaks at M ~ 2 x 10 14 Mq, i.e., for haloes of Lagrangian radius 
r ~ 9/i - 1 Mpc. In Fourier space this corresponds to a wavenumber k ~ ir/r = 0.35/iMpc -1 . If 
we assume that this stochastic term induces a Gaussian smoothing in the resulting Eulerian 
space density field, we conclude that it produces a 1 % change in the density power spectrum 
at k ~ 0.2/iMpc _1 and will thus contribute to the Eulerian EFT sound speed c s . The 
time dependence can be introduced into the halo model by scaling the variance as a(M) —> 
D(a)cr(M). The resulting u,r toch h then scales as D i e . For this scaling we employed a Sheth- 
Tormen mass function [27]. 

Fitting function The stochastic term can be accounted for by a simple fitting function 
that is constructed from the following requirements: 


- 31 









• on large scales it scales as fc 4 ; 


• on small scales it scales as k m with m < — 1; 


• it integrates to the stochastic displacement dispersion cr 2 t0( , h . 

One function that satisfies these constraints and smoothly transitions between the low- and 
high-A; regimes is given by 


k 4 P 


f k\ A 2 1 2 4 t r 3 / 2 r(n) 

W astoch (l + (k/k s y^ k s r(n — 5/2) ’ 


(5.10) 


where n > 2.5. As we show in Fig. 23 this functional form provides a very good description of 
the measured stochastic term at all redshifts once the scalings cr s toch oc D 15 and k s oc D~ 15 
are accounted for. For the scale we employed k s ~ k n \ = 0.26hMpc _1 . 


5.2 Origin of the stochastic term 

What is the origin of this observed stochastic term </> s toch? Such a term is expected in the 
framework of the EFT, due to the the small-scale fluctuations which are not amenable to 
perturbation theory and are treated as a random noise. This term could also arise from 
higher order LPT terms which we haven’t computed. Finally, it could also be a mere artefact 
of numerical errors in the A-body simulation. 

In order to exclude the possibility that this stochastic term is simply due to numerical 
errors, we have run simulations with two different box sizes. The stochastic terms .P s toch iden¬ 
tified at redshift z = 0 in these two simulations agree, as can be seen in Fig. 23. In this context 
we should mention that the change in the simulation box size also changes the wavenumber 
at which we have to cut our initial power spectrum in order to avoid aliasing. Thus, while 
the L simulation has fc max = 0.6 h^ 1 Mpc the M simulation has fc max = 2.4 h^ 1 Mpc. This 
difference in the cutoff wavenumber translates into slightly different LPT contributions and 
slightly different transfer functions. Yet, at the scales where we trust the simulations, the 
power spectrum of the stochastic term does not change. We interpret this observation as a 
strong indication of this term being truly related to highly non-linear motions that are uncor¬ 
related with the perturbation theory prediction. Virialized motions are an example of such 
motions uncorrelated with the LPT terms, since LPT is not able to capture shell crossing. 

We also considered transients in the IV-body simulation [23] as a possible source of 
stochastic term. Indeed, the A-body simulations are initialized with 2LPT (at redshift z m a = 
99), so we expect the simulations to contain transient errors in the third and higher order 
displacement fields. This error can be estimated in LPT in terms of an error on the growth 
factor of each LPT term (see App. D): fa and fa are not affected, fa is underestimated by 
2 x 10 , the £ 2,2 and Ad 1 ’ 1,2 contributions to </>4 are underestimated by 10~ 6 while the £ 3,1 

contribution to fa is underestimated by 4 X 10 -4 . These errors are too small to account 
for the measured stochastic term, and they are clearly correlated with the LPT terms, so 
they should be automatically corrected for by the transfer functions in tLPT. In conclusion, 
transients cannot account for the stochastic term. 


5.3 Link with the curl part of the displacement 


Besides the divergence, we can also consider the curl of the displacement field, 
correlators we have 


((V x t/>)i|(V x ip)i) = (AwilAwi) , 


For the 
(5.11) 
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Figure 25. Power spectra of the stochastic displacement divergence V • i/’stoch (red dashed), the full 
curl V x i/i (black solid) and the LPT prediction for the curl (gray). From top to bottom we show 
measurements and predictions for z = 0, 2 = 1 and 2 = 2. At redshift 2 = 0 the measurement of the 
curl exceeds the LPT prediction by roughly one order of magnitude but agrees with the divergence of 
the stochastic displacement on large scales. The dot-dashed lines at 2 = 1 and 2 = 2 show the sum of 
the LPT curl and the stochastic part of the displacement divergence. 


where we used Coulomb gauge. Thus, we have in Fourier space 

-Pvxi/’ = ((& x ij))i\(k x 1 / 7 );) = k 4 {oji\uii) . (5.12) 

Correlating the curl components of the field itself, we have 

(ipc,i IV’c.i) = ((V x u>)j|(V x a ;)i) (5.13) 


or in Fourier space 

P^ c = ((^ ^ ^ ^)z) = = 

=k 2 (ui\u>i) , (5-14) 

where we again used the Coulomb gauge condition kiUJi = 0. The first quantity is easier for 
measurements in simulations while the latter is more closely related to perturbation theory. 
As we saw above, their power spectra are related by multiplicative factors of k. 

We measure the curl power spectrum T’vxvE' hr both the M and L simulations and find very 
good agreement where the spectra overlap. In agreement with [28] we find that the curl 
power spectra scale as D 10 . As shown in Fig. 25, at redshift 2 = 0 the measured curl power 
spectrum exceeds the LPT prediction from P 3 C 3 C by about a factor of 10 on large scales. 
This indicates that the curl component of the displacement is completely dominated by small 
scale, non-perturbative virialized structures. Since the perturbative prediction scales as D 6 , 
we expect the perturbative and non-linear curl power spectra to be of the same order only at 
2 = 2. Indeed we see a slight upturn of the curl power spectrum at low k at 2 = 2, where the 
measured curve approaches the perturbative prediction. 

Since the leading curl contribution is fully non-perturbative, it is interesting to compare its 
power spectrum to the stochastic term of the displacement divergence, which we argued to 
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be sourced by virialized motions. Amazingly, both the large scale amplitude and the redshift 
scaling of the curl are in very good agreement with the scalar stochastic term. This further 
solidifies our confidence in having identified the true stochastic term arising from virialized, 
non-perturbative motions. 

6 Conclusions 

In this paper we studied Lagrangian Perturbation theory and its regularization in the effective 
held theory approach. We numerically evaluated the LPT displacement fields up to fifth order 
on a regular grid that has the same phases as a corresponding suite of A-body simulations. 
This allowed us to test the 1-loop LPT and EFT at the level of the displacement held itself, 
which is a much more stringent test of the perturbative approach than matching the power 
spectrum only and also reduces the cosmic variance significantly. As an added beneht, we 
were able to evaluate LPT terms (such as P 55 , which is formally a 4-loop term) that would 
be barely tractable with a Fourier space loop calculation. 

We verified that the LPT expansion is well-behaved up to third order, in that the error on 
the displacement held decreases from 1LPT to 2LPT, to 3LPT, on scales k < O.lhMpc -1 . In 
doing so, we highlighted the importance of the second order displacement cj) 2 ? which is crucial 
in reducing the error power spectrum (even more so than (^ 3 ), despite its small contribution 
to the non-linear power spectrum (P 22 -P 13 ). We found that 1-loop LPT provides a 1%- 

accurate description of the non-linear power spectrum up to k = 0.05/iMpc -1 . 

We carefully tested the 1-loop EFT. We validated the simple scalings of the various LPT 
terms and the EFT counterterm up to 1-loop, estimated for an EdS Universe. We reliably 
detected the leading EFT counterterm ak 2 (j) 1 with a non-zero coefficient a on large scales 
(k < 0.03/rMpc"” 1 ). The EFT provides a significant improvement over LPT, with a better 
convergence up to 1 -loop, and a threefold increase in the maximum wave vector where the 
non-linear power spectrum is reproduced to 1%: the 1-loop EFT power spectrum is accurate 
to 1% up to k = 0.15/iMpc -1 . The EFT also reduces the error power spectrum by up to a 
factor of two around k = 0.1/iMpc -1 . The measured time-dependence of the EFT coefficient 
a matches reasonably well the one expected from simple scalings in an EdS Universe. We 
showed that the EFT is close to optimal, in the sense that it achieves similarly small error 
power spectrum P error and similarly accurate non-linear power spectrum as the 3tLPT and 
higher order tLPT models do. 

By looking at the displacement held and not only its power spectrum, we were able to 
detect a stochastic term </> s toch that is uncorrelated with the LPT terms and seems to constitute 
a limit to the reach of LPT and EFT. Its power spectrum is slightly smaller than that of q 
but larger than that of and the corresponding root-mean-square displacement in real 

space is consistent with the typical random motion within halos. Its scale-dependence on 
large scales matches what is expected of a momentum-conserving displacement, consistent 
with the stochastic term expected in the EFT framework. However, the size of this effect 
seems larger than the naive expectation based on power law Universe scalings. To the best of 
our understanding this term is not due to numerical errors or transients in our simulations. 
For example, its amplitude matches in simulations that have very different resolutions. If not 
modeled explicitly, this term would imply a limit to the accuracy to which the displacement 
held can be described by LPT or EFT, corresponding to a 1% error on the non-linear power 
spectrum at k = 0.2hMpc -1 at 2 = 0. In this case, there is no beneht in going to higher 
order than third order in the displacement with the corresponding one-loop counterterm. 
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In particular, no tLPT-like model fit to the displacement can predict the non-linear power 
spectrum to better than 1% beyond k = 0.2 h/Mpc, unless it overfits the power spectrum. 
This conclusion extends the result from Fig. 20 to all orders in tLPT, and not only at one- 
loop. However, one might be able to model the power spectrum of this stochastic term with a 
few free parameters or a fitting function and thereby extend the range of validity of the EFT 
displacement power spectrum. 

For practical calculations of the density held, Eulerian or Standard Perturbation Theory 
and the corresponding Effective Field Theory are very useful. Performing a similar test at 
the level of the fields rather than power spectra would be desirable, but is complicated by the 
decorrelation due to long wavelength motions. Thus the long wavelength displacements need 
to be resummed [29]. We will address this issue in a forthcoming paper. 

Besides testing the EFT approach, we also pushed the simulations to their limits. The 
EFT, or more generally its underlying symmetry arguments, predict how corrections to the 
linear power spectrum should behave on large scales. These constraints are not fully satisfied 
by the simulation measurements at the sub-percent level. In this context, we pointed out 
that the measurement of the a coefficient on large scales is very difficult, because it is highly 
sensitive to systematic errors in the simulations. 
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A Numerical tests 

We estimate the systematic error inherent to our simulations by varying several parameters 
of the Gadget IV-body simulation and study the deviations between these simulations and 
our fiducial case. The corresponding systematic error on the power spectrum is shown in 
Fig. 26. In particular we estimate the cross power and auto power of the differences between 
two simulations run with parameter choices A and B 

Pa-BxB = ({&A —<Pb)\4>b) , Pa-b = ((4>a ~ 4>b) 2 ) ■ (A.l) 

It turns out that the results are very sensitive to the number of grid cells that Gadget uses 
to calculate the long range potential. This parameter is called PMGRID and has to be set in 
the Gadget makefile. We had initially set AVm.f = N c , which we shall refer to as the fiducial 
case. We also changed this parameter to Apm,pm = 2A C or IVpmpml = 3/4A c . As was 
noted previously in [30], the PMGRID parameter has a big and non-monotonic effect on the 
matter density power spectrum. As shown in Fig. 26, we confirm a similar behavior for the 
displacement divergence. In particular we find that the auto power of the difference between 
the PM and PML cases is smaller than the respective errors between the F and PM/PML 
cases. We also considered a simulation with an improved time stepping and error tolerance 
(denoted HR in the Figure) and find a considerably smaller impact on the displacement power 
spectrum. We also consider the time dependence of the error power spectra and find that 
both the auto and cross power spectra scale roughly as D 2 . We find fitting functions that 
roughly reproduce the shape and amplitude of the errors. These fitting functions are used 
in the main text to assess the systematic error on the measurements of the EFT coefficients 
and the stochastic term. In particular, the cross power of the error and the field Pa-BxB 
will affect the transfer functions and the EFT coefficients, while the auto power Pa-BxA-B 
affects the stochastic term. 

B Perturbative Displacement Fields 

In this Appendix we will rederive the equations for the scalar and vector components of the 
displacement, as they will be needed to calculate the higher order solutions in presence of 
EFT counterterms in App. E. 

Scalar component of the displacement: For the scalar equation we take the Eulerian 
derivative of the EoM, yielding 

r\ q 

— + nty = —Aip = --n 2 n m 6 (b.i) 
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Figure 26. Power spectrum of the difference displacement divergence between simulations with 
variant parameters. The solid points show the auto power of the error, while the crosses show the 
cross correlation between the error and one of the displacement divergencies. Right panel: Time 
dependence of the error from measurements at z = 0 and z = 2. Both errors on the auto and the 
cross power scale as D 2 . 


From the mapping between Lagrangian and Eulerian space x = q + ip, we have that 5 = 
1/J— 1, where J is the determinant of the Jacobian matrix 

Aij = || = C + (B- 2 ) 

If not explicitly mentioned otherwise, partial derivatives are with respect to the Lagrangian 
coordinate q. The determinant can then be related to the invariants of the displacement field 

J — Det [A] — —\£iji£stu-AisAjtAi u 

— 1 + 1pi,i T ~ + T^^ijRstu^Pi.sPj.t^Pl.u (B.3) 
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We will need the mapping from Eulerian to Lagrangian derivatives 
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For the EoM we now have 

JAjl {Ppi,j + = 2^m^ 2 (/C + C + AJ) 

Let us define the time derivative operator 
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We can further simplify this equation by going to an EdS Universe and keeping only the n-th 
order displacement held on the left hand side 
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As a last step we can symmetrize the right hand side and collect the invariants C and M. 
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In the second part of the above equation, we have restored the full index structure, which 
will become useful for writing the recursion relations in Fourier space later on. 


Vector component of the displacement: For the curl part of the displacement held, we 
start again from the equation of motion: 


Xi + T-Lxi 
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dx{ 
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The Eulerian gradient on the r.h.s. can be converted to a Lagrangian gradient by multiplying 
the equation by the Jacobian matrix Ajj = ^: 


dxi 

dqj 


(xi + T-Lxi) 
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With integration by parts with respect to time and space, this equation can be reexpressed 
as: 



d Qj 1 d qj \ 2 )■ 


(B.14) 
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The point of this transformation is that the r.h.s is now a pure Lagrangian gradient, and 
therefore its Lagrangian curl vanishes: 


^Imj 


d ( d 


dq m \dr 
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(B.16) 


dqj dq m 

This differential equation in time is readily integrated, since the initial condition vanishes: 

dxi dii 


___ - i' -^ 

Qmjv; t; — 

oqj dq m 

which can be rewritten in vector notation as: 

^qXi X V qXi — 0. 

Substituting x = q +ip then leads to: 

V q X1p = V q'lpi X Vqlpi. 


(B.17) 
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Assuming that the time-dependence of the displacement follows that for an EdS universe, 
this yields the configuration space recursion relation for the LPT curl part: 


71—1 
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Recursion relations in Fourier space: Finally, we deduce the Fourier space recursion 
relations for the kernels of the scalar l n and vector part t n of the displacement. Both are 
sourced by the total displacement S n = ikL n + ik x T n , which leads to a coupling between 
scalar and vector modes starting from third order. 

The scalar kernels L n , n > 1 are readily obtained by Fourier transforming Eq. (B.ll): 
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with starting conditions L\ = 1/k 2 and T\ = T 2 = 0. Here, Det [a, 6, c] = a ■ (b X c) stands 
for the determinant of the matrix, whose columns are the vectors a, b and c and we defined 
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pI™ = pi + ... + p m - For the vector kernel T n , we get from Eq. (B.20) 
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Similar relations were recently given by [31] . For the purely scalar part Eq. B.21 simplifies 


to 3 4 
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This equation can be used to calculate the full scalar displacement up to third order and the 
dominant part of the fourth order displacement. 

In order to validate the LPT code on the simulation grid, we compare it to the analytical 
1-loop and 2-loop predictions. They can be computed using the following explicit relations 


3 These expressions correspond to Eqs. (67-69) in [31] once the following mapping has been used 
In.Matsubara = S n , S n , Matsubara = ik 2 L n and T„, M atsubara = ik 2 T n , where k = p|? is the sum of the mo¬ 
menta in the kernel. 

4 Here we are using that 
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corresponding to the diagrams in Fig. 1 
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We find a good agreement between onr code on the grid and the loop calculations, as 
seen in Fig. 27. 



Figure 27. Comparison between the LPT on the grid and the analytical one, for the linear power 
spectrum, 1-loop terms, and some 2-loop terms. This is a test for our expressions of 4>^ l \ <\P'\ 
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For definiteness, we have for the one and two loop LPT contributions to the power 
spectrum 

Fuoop = 2 Pia + P 22 , (B.25) 

-P2-I00P = 2P15 + 2P 2 4 + P33-1 + P33-11 ■ 

The calculation on the grid has effective cutoffs: the fundamental wave vector fcf = 
(27r)/Lbox is the minimum non-zero wave vector present, and if no explicit upper cutoff is used, 
the Nyquist wave vector k^ y is the highest wave vector present. In order to get agreement 
between the theory LPT and the calculation on the grid, one has to impose these same cutoffs 
to the loop integrals. Since the LPT expansion is a non-linear calculation, aliasing is also 
an issue, as explained in the main text. Finally, reaching an agreement at the percent level 
between theory LPT and LPT on the grid requires binning the theory power spectra with the 
same bins as for the power spectrum estimation on the grid. 


C Cosmic variance 


In this appendix we estimate the uncertainties on our measurements due to cosmic variance, 
in order to show the difficulty of measuring the EFT coefficient at the level of the power 
spectrum (instead of the field itself), and to understand the scatter in our measurements of 
the cross power spectra between the displacement from simulation and from perturbation 
theory. 

The variance of our measurements can be estimated from their scatter across the various 
realizations. On the other hand, these variances can be predicted as: 

Cov (Pah, P cd ) = — - [PacPbd + PadPbc] + TjTabcdr (C.l) 

™ modes v 


where N modes is the number of modes in each fc-bin, V is the volume of the box, and P a b cd is 
the trispectrum. The first term is the Gaussian covariance, while the second term is the non- 
Gaussian term, which would vanish if the fields of interest were Gaussian. In what follows, 
we only evaluate the Gaussian covariance, although we will keep the trispectrum terms in the 
equations for completeness. 

In particular: 
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where r a b is the correlation coefficient r a b = 
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^ 1. A simple consequence of Eq. (C.2) 
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spectrum on large scales, and it corresponds to a relative error of ~ 10% on the largest scales 
of our simulations (k ~ 5 x 10“ 3 /iMpc -1 ), as shown in Fig. 28. Such a high uncertainty would 
make it extremely difficult to measure the EFT coefficient a from matching the simulation 
power spectrum to a theory power spectrum on large scales. 

This uncertainty can be significantly reduced by matching ratios of power spectra from 
the same simulation to theory predictions. This can be understood as follows. The uncertainty 
on a ratio of power spectra can be computed from Eq. (C.l) as: 
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Figure 28. Relative scatter of the various power spectra and ratios of power spectra, due to sample 
variance. On large scales, considering ratios of power spectra cancels the sample variance by up to a 
factor of a hundred. 



Figure 29. Relative uncertainty on P n u due to cosmic variance (Gaussian term only), showing that 
the relative uncertainty is much greater for i even. This explains the scatter seen in the measurements 
of Fig. 16. 


Fig. 28 shows that considering ratios of power spectra reduces the cosmic variance from 10% 
to 0.05% at k ~ 5 x 10~ 3 /iMpc -1 . Such an uncertainty would be marginally sufficient to 
allow our measurement of the EFT coefficient a. 

Another consequence of Eq. (C.2) is that the relative uncertainty on P n n is much greater 
when i is even, since r n p is much lower for even values of i. Indeed, r n p = —^=== ~ rp for i 


odd, while r n u — y T^ r 2i ^2 i for * even. This is shown in Fig. 29, and explains the scatter 

in the measurements seen in Fig. 16. 


D Transients 

In this Appendix, we consider the imprint of a finite starting redshift of the simulation on 
its late times results. These so called transients were studied in perturbation theory [32] and 
simulations [23] at the level of 1LPT initial conditions. We will review the derivation and 
extend it to higher orders. 
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The equation of motion can be written in vector notation for \ = (ipi i, rji^) = (^i i,dipi i/d In a) 
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Here we already brought the IC^d = term arising from the determinant Eq. (B.3) to 
the left hand side. The source term is thus composed of the C^ and terms in the 

determinant and the derivatives of the lower order terms from the left hand side. 

The above equations can be summarized in vector notation as 


with the solution (see e.g. [32]) 
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Let us assume that we use Zel’dovich initial conditions, i.e., the 1LPT part is correct at all 
times, but the second and higher order solutions are zero until some initial time dj. Then the 
second order source term is given by 
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The first term in the brackets gives the fastest growing solution, whereas the second and third 
terms lead to corrections that decay as with increasing a. For a starting redshift z\ = 99, 
i.e. a; = 0.01 this corresponds to a 1.4% correction to the second order displacement field at 
redshift z = 0. 

Let us now consider the case where we use 2LPT to set up initial conditions at a;, i.e., 
first and second order displacement fields are correct at all times and there are transient effects 
on the third and higher orders. This is the case that is relevant for our suite of simulations. 
The third order source term is correctly given by the growing modes 

£(3)(a) = a 3 [—3M 1 ’ 1 ’ 1 - IOC 1 ’ 2 ] , (D.9) 
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/o') 

which using © (a < a;) 


0 readily leads to the following solution for a > a; 



(D.10) 


Here the first term in brackets gives the fastest growing mode and the corrections amount to 
2 x 10 -4 at 2 ; = 0 for z\ = 99. 

The fourth order source term has parts that are correctly predicted by LPT, but also decaying 
mode corrections that lead to corrections in the source term. It can be readily obtained by 
using the solution for the third order held including transients Eq. (D.10) in Eq. (B.9) 
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Thus we have for the fourth order solution including transients 


V’i? =« 4 
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The relative error for C 2,2 and M 1,1,2 is 2 X 1CP 6 and the relative error on £ 1,3 is 1CP 4 at 
2 ; = 0 for Zi = 99. The enhanced error on £ 1,3 is due to the presence of transients in the third 
order contributions to the source term. This term decays slower than the transients on the 
other components of the fourth order displacement held by one power of the expansion factor 
a. 


E Solution with Source Terms 


Let us explicitly consider the time dependence of the EFT counterterms that arises from a 
given time dependence of the source term on the right hand side of the equations of motion 
(EoM). For convenience, we factor out 3/2Q m 'H 2 from the source terms, such that they 
appear similar to the gravitational potential and such that that the EFT parameter has units 
of length 2 even at the EoM level 

4> {i \k) + n^ l \k) - ^n m n 2 ^\k) = ^n m n 2 pk 2 ^ 1 \k ). (e.i) 


From now on we will work in Einstein-de-Sitter, where fl m = 1 and postulate a power law 
time dependence for the counterterms in the equation of motion, for instance = /©a 7 . It is 
now convenient to rewrite the equations using In a as the time parameter: 


a) = ^/3q a 7 ©^ 1 -^© a) with V 
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Note that here we redefined the operator T>, to also contain the full homogeneous part. The 
Green’s function associated with V is given by 
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This yields for the integration of a power law source 


[ din a G(a,d)a w = - - ^ - -a w (E.4) 

Jo K 1 (w-l)(3 + 2 w) V ' 

and hnally allows us to relate the EFT coefficient a of the solution to the one of the equation 
of motion 

4>^\k,a) = a(a)k 2 ((M\k,a) with a(a ) = —-r/?oo 7 - (E.5) 

7(5 + 27) 

Let us now solve the equation at next-to-leading order in the counterterms, but neglecting 
the explicit second order source terms 772,?: for the time being. 

V(jJ 2 Xk,a) = ^/? 0 a 7 /cV 2) (fc,a) + p ^(7 + l f + (E.6) 


The first term arises from considering the second order LPT held in the leading EFT coun¬ 
terterm and the second source arises from coupling the linear held and the leading order EFT 
solution using the the quadratic kernel. More explicitly we can write 
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Integrating the source with the Greens function, we hnally have 
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We thus see that expressions in the form of the 772,2 and 772,3 counterterms are generated by 
the time integration of the leading order counterterm, but not 772,i- Part of the counterterm 
<; 1 u 2 ' appears in 77,,, to cancel the UV mistake in P 24 . For this part of the counterterm, the 
time-dependences of 77,3 and P 24 have to match match, which implies 7 = 2. 


47 



